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ABSTRACT 

Recent discoveries of circumbinary planets by Kepler mission provide motivation for understanding 
their birthplaces — protoplanetary disks around stellar binaries with separations < 1 AU. We explore 
properties and evolution of such circumbinary disks focusing on modification of their structure caused 
by tidal coupling to the binary. We develop a set of analytical scaling relations describing viscous 
evolution of the disk properties, which are verified and calibrated using ID numerical calculations with 
realistic inputs. Injection of angular momentum by the central binary suppresses mass accretion onto 
the binary and causes radial distribution of the viscous angular momentum flux Fj to be different from 
that in a standard accretion disk around a single star with no torque at the center. Disks with no mass 
accretion at the center develop Fj profile which is flat in radius. Radial profiles of temperature and 
surface density are also quite different from those in disks around single stars. Damping of the density 
waves driven by the binary and viscous dissipation dominate heating of the inner disk (within 1-2 
AU), pushing the iceline beyond 3-5 AU, depending on disk mass and age. Irradiation by the binary 
governs disk thermodynamics beyond ^ 10 AU. However, self-shadowing by the hot inner disk may 
render central illumination irrelevant out to 20 AU. Spectral energy distribution of a circumbinary 
disk exhibits a distinctive bump around 10/im, which may facilitate identification of such disks around 
unresolved binaries. Efficient tidal coupling to the disk drives orbital inspiral of the binary and may 
cause low-mass and relatively compact binaries to merge into a single star within the disk lifetime. 

We generally find that circumbinary disks present favorable sites for planet formation (despite their 
wider zone of volatile depletion), in agreement with the statistics of Kepler circumbinary planets. 

Subject headings: planets and satellites: formation — protoplanetary disks — stars: planetary systems 


1. INTRODUCTION. 

Theory of planet formation has experienced a wave of 
recent activity motivated by the large increase in the 
number of known extrasolar planetary systems. Some 
of these systems are rather unusual, which is nicely il¬ 
lustrated by the recent discoveries of circumbinary plan¬ 
ets (affectionately termed ’’Tatooines”) around about a 
dozen of eclipsing binaries by the Kepler mission (Doyle 
et al. 2011; Welsh et al. 2012; Orosz et al. 2012a; 
Schwamb et al. 2013). Central binaries in these systems 
are composed of main sequence stars with masses around 
or below Mq. They have semi-major axes ab = 0.08 — 0.2 
AU and some of them are quite eccentric (binary Kepler- 
34 has eccentriity Cb = 0.52). Characteristics of plan¬ 
etary orbits in these systems are measured with good 
accuracy thanks to the exquisite timing precision of the 
Kepler satellite, and they are found to be coplanar with 
their binaries to within a few degrees. 

Further, presence of planets orbiting eclipsing post¬ 
common envelope binaries has been suspected based on 
eclipse timing of these systems. A noteworthy example 
is the evolved binary NN Serpentis (Beuermann et al. 
2010; Marsh et al. 2014) for which it has been argued 
that if real, the planets would have formed from the ma¬ 
terial ejected during the common envelope phase (Mustill 
et al. 2013; Volschow et al. 2014). 

These discoveries show that planet formation in cir- 
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cumbinary disks is feasible, despite the complications 
related to the gravitational perturbations induced by 
the central binary. Recent theoretical efforts have ad¬ 
dressed some of these issues (Paardekooper et al. 2012; 
Alexander 2012; Meschiari 2012; Gong et al. 2013; Fou- 
cart & Lai 2013; Pelupessy & Portegies Zwart 2013; 
Rafikov 2013a; Clanton 2013; Martin et al. 2013; Silsbee 
& Rafikov 2015). They demonstrated, in particular, that 
an instrumental part of the theory behind the circumbi¬ 
nary planet genesis is the understanding of the prop¬ 
erties and evolution of the circumbinary protoplanetary 
gaseous disks. Disk properties determine the dynamics 
of planetary building blocks — planetesimals — in many 
ways. They set not only the efficiency of gas drag but 
also the precession rates of both planetesimal orbits and 
the binary itself (Rafikov 2013a; Silsbee & Rafikov 2015). 
Even for post-common envelope binaries, it has been con¬ 
tended that a significant fraction of the ejected material 
forms a circumbinary disk (Kashi & Soker 2011). 

Circumbinary disks provide a good testing arena for 
ideas about planet formation, in part because the known 
systems with circumbinary planets are well constrained, 
and also because the presence of the binary has impor¬ 
tant consequences for the evolution of the disk. The bi¬ 
nary torques are expected to clear out the inner disk re¬ 
gion (Artymowicz & Lubow 1994; Pelupessy & Portegies 
Zwart 2013), significantly reducing the mass supply to 
the binary (MacFadyen & Milosavljevic (2008), but see 
D’Orazio et al. (2013)). Moreover, the transfer of angu¬ 
lar momentum from the binary to the disk changes the 
structure and dynamical evolution of the disk on large 
scales. Understanding of the interplay between these 
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complicated processes is necessary to gain insight into 
the circumbinary planet formation. 

The goal of this study is to explore general proper¬ 
ties and and main evolutionary features of circumbinary 
protoplanetary disks around stellar binaries. We focus 
on understanding the differences between the character¬ 
istics of the circumbinary disks and their more stan¬ 
dard counterparts around single stars. In doing this 
we fully account for the disk thermodynamics influenced 
both by viscous heating and irradiation by the central 
binary as well as by tidal interaction between the disk 
and the binary, through the dissipation of the density 
waves launched by the binary in the disk. In the course 
of our study we provide a semi-analytical description for 
the evolution of disk properties that we test and cali¬ 
brate using detailed numerical calculations. Our results 
are then used to assess the implications for circumbinary 
planet formation, and for the evolution of the central 
binary itself. 

This paper is organized as follows. After describing our 
general setup in §2, we cover the basic of the viscous evo¬ 
lution of circumbinary disks in §3. In §4 we describe our 
treatment of disk thermodynamics. In §5 we derive a set 
of analytical scaling relations describing the viscous evo¬ 
lution of the disk properties. We then verify these results 
numerically in §6 (our numerical approach is outlined in 
Appendix A). Spectral differences between the circumbi¬ 
nary and circumstellar disks are described in §7. In §8 
we discuss the role of different heating sources (§8.1), 
effect of accretion onto the binary (§8.2), orbital evolu¬ 
tion of the binary due to the tidal coupling with the disk 
(§8.3), dead zone (§8.4) and iceline (§8.5), limitations of 
our models (§8.6) and provide comparison with the ex¬ 
isting work on circumbinary disks (§8.7). We cover the 
implications for circumbinary planet formation in §9 and 
conclude with a brief statement of our findings in §10. 

2. GENERAL SETUP. 

We focus on a particular case of a thin disk (we 
disregard its vertical dimension) coplanar with the bi¬ 
nary. The binary has semi-major axis ab, total mass 
Me = Mp + Ms {Mp and Mg are the masses of the pri¬ 
mary and secondary); the mass ratio of its components 
IS q = Mg/Mp < 1. 

Our description of the circumbinary disk assumes it 
to be axisymmetric, with all characteristics — surface 
density E, midplane temperature T, etc. — to be func¬ 
tions of radius r only. Thus, we neglect the short-term 
variability of the disk properties caused by the orbital 
motion of the binary and focus only on the long-term, 
time-averaged effect of the binary on the disk. This in¬ 
fluence comes mainly in two flavors. 

First, the gravitational potential in which the disk or¬ 
bits is time-variable, which drives density waves in the 
disk. Although we do not explore the two-dimensional 
structure of these perturbations, we do account for this 
tidal coupling by including the associated angular mo¬ 
mentum and energy injection in the disk in our calcu¬ 
lations (§3 and §4.1). But to zeroth order, we model 
the disk as orbiting in the potential of a point mass M^ 
centered on the barycenter of the binary. 

Second, illumination of the disk by the binary plays im¬ 
portant role in setting its thermal structure (§4.1, 8.1). 
Orbital motion of the binary leads to periodic variations 


of the stellar flux impinging on disk surface, the effect of 
which has been explored in Clanton (2013) and Bodman 
& Quillen (2015). We neglect this variability of irradia¬ 
tion and study only its time-averaged effect in this work. 


3. VISCOUS EVOLUTION OF A CIRCUMBINARY DISK. 

Evolution of the circumbinary disk is driven by both 
viscous stresses and angular momentum injection by the 
binary. It is well known that one-dimensional (only in 
r) viscous evolution of the disk surface density E with 
external sources of angular momentum can be described 
by a single diffusion equation (Papaloizou & Lin 1995): 
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Here H is the angular speed and I = il{r)r^ 

is the specific angular momentum in the central binary 
potential. Following standard approach we model viscos¬ 
ity V via the a-ansatz (Shakura & Sunyaev 1973) 
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where T and Cg are the disk (midplane) temperature and 
sound speed, and /i is the mean molecular weight of the 
gas. 

The injection of the angular momentum by the binary 
into the disk can be characterized via the specific angu¬ 
lar momentum injection rate A(r), with dimensions [cm^ 
s“^] in CCS units. It is defined as the amount of angu¬ 
lar momentum transferred by the binary to the disk at 
radius r per unit time t and unit disk mass. 

In this work we follow Armitage & Natarajan (2002) 
and adopt 


A(r) = sgn(r - Ob)/ 
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where f is a dimensionless parameter that is chosen as 
described in Appendix A. This prescription, namely the 
power law scaling with \r — ab\, is motivated by the an¬ 
alytical calculations of Lin & Papaloizou (1979b) and 
Goldreich & Tremaine (1980). Recent numerical work 
(MacFadyen & Milosavljevic 2008; Cuadra et al. 2009) 
shows that torque density in a disk around a high mass 
ratio binary exhibits considerably more complicated, os¬ 
cillatory pattern at the cavity edge, rather than following 
a simple scaling (3). The analytical work of Rahkov & 
Petrovich (2012); Petrovich & Rahkov (2012) in the low-g 
limit also suggests that near the gap edge the prescrip¬ 
tion (3) may be incomplete at best. Non-local nature of 
the density wave damping (Goodman & Rahkov 2001; 
Rahkov 2002) is also not accounted for by equation (3). 
Nevertheless, because of its simplicity and also for rea¬ 
sons outlined in Appendix A we use the prescription (3) 
in this work. It should also be remembered that the last 
term proportional to A should be absent in equation (1) 
in the case of a circumstellar disk. 

Torque exerted by the binary on the disk clears out 
a cavity in the disk center provided that the mass ratio 
of the binary components is not too small (Artymowicz 
& Lubow 1994). For q ^ (0.1 — 1) the cavity radius 
Tc is about twice the binary separation (MacFadyen & 
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Milosavljevic 2008), although also shows some depen¬ 
dence on the binary eccentricity (Pelupessy & Portegies 
Zwart 2013). At the cavity edge surface density falls off 
quite abruptly within a narrow interval of radii. This is 
caused by a very steep decay of the torque density with 
r (Lin & Papaloizou 1979b; Goldreich & Tremaine 1980; 
Petrovich & Rafikov 2012) — a point that is further dis¬ 
cussed in Appendix A. As a result, small reduction in 
r leads to large increase of A preventing viscous inflow 
of material from the disk. However, it is worth point¬ 
ing out that non-axisymmetric gas motions at the cavity 
edge (not captured by the ID approach based on equa¬ 
tion (1)) do give rise to some accretion from the disk 
into the cavity, albeit at a significantly reduced rate 
than in the case without a central binary (MacFadyen & 
Milosavljevic 2008). 

Rafikov (2013c) has shown that properties and evolu¬ 
tion of circumbinary disks can in many ways be better 
understood if instead of E one uses the viscous angular 
momentum flux Fj defined as 

Fj = —27ri^Er^-— = SttvEI, (4) 


where the last equality is for a Keplerian disk with 
Q oc . By definition Fj is the total viscous torque 
exerted by the inner disk on the outer disk at a given 
radius. 

Additionally switching from r to the specific angu¬ 
lar momentum f the evolution equation (1) can be re¬ 
written as 
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is the diffusion coefficient. Note that Dj depends on Fj 
if 1 / depends on E. 

In terms of the new variables I and Tj, the mass ac¬ 
cretion rate (M > 0 for mass inflow) is given by 
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This formula shows that the binary torque suppresses 
mass accretion into the cavity as A > 0 for r > ah- 


3.1. Disk evolution away from the cavity edge. 

The aforementioned steep dependence of A on r im¬ 
plies that the binary torque becomes insignificant for the 
disk evolution outside some radius r\, which is not very 
different from the cavity radius Vc — radius where E 
(and Fj) becomes small. Thus, it is convenient to sep¬ 
arate the disk into two parts, as illustrated in Figure 1. 
One is the edge region {vc ^ r < r\) in which A is large 
and where essentially all of the angular momentum and 
energy carried by the binary-induced density waves get 
dissipated. Second is the bulk of the disk in which A is 
negligible and the surface density (and Fj) evolution is 
determined solely by viscous stresses. 

The utility of switching to Fj and I variables can be 
best appreciated in the bulk of the disk, r > r\, outside 
the cavity edge region. As A —0 there, the equation (5) 


Cavity Disk Edge Bulk of the Disk 



Fig. 1. — Schematic illustration of the behavior of the disk surface 
density E, angular momentum flux Fj and specific torque density 
A in different parts of the circumbinary disk. Vertical dashed lines 
show the binary semi-major axis radius of the inner cavity rc 
and the radius beyond which the binary torque density can 
be neglected. Key disk regions discussed in the text — cavity 
(t’ ^ Tc), disk edge (vc ^ r < ^a), and bulk of the disk (r > r\) 
— are indicated. The case of non-accreting binary (M^ 0) is 

illustrated, so that Fj is radially constant in the bulk of the disk. 

transforms to a particularly simple form (Lynden-Bell & 
Pringle 1974; Filipov 1984; Lyubarskij & Shakura 1987; 
Rafikov 2013c) 


d fFj\^ d^Fj 
dt\Dj) dP 


( 8 ) 


Also, M {l,t) = dFj/dl in this case, according to equa¬ 
tion (7). 

Equation (8) must be supplemented by a boundary 
condition (BC) at the inner boundary of the main disk 
region, i.e. at r = ta, which characterizes the suppres¬ 
sion of the mass inflow by the binary torque. We express 
this condition via the constraint on the mass accretion 
rate of gas from the disk into the cavity Mf,. We assume 
Mb to be set by the processes happening in the edge re¬ 
gion, at r < ta (MacFadyen & Milosavljevic 2008). Its 
value is set by the presice form of A(r) and can be non¬ 
zero in general. 

Given the small radial extent of the edge region, it 
contains rather small amount of mass and one can safely 
assume that M at the inner boundary of the bulk of the 
disk, at rA, is just equal to Mb. This means that a proper 
inner BC for equation (8) is 


dFj 


dl 


K^a) 


Mb. 


( 9 ) 


When studying global disk evolution at r ^ ta, one can 
simply assume this inner BC to be imposed at r —>■ 0 
9 ). 

In this work we will mainly work with the no-inflow 
BC, such that Mb = 0. 
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3.2. (Quasi-)steady state cireumbinary disks. 

We now use the framework outlined above to under¬ 
stand the steady state structure to which a cireumbinary 
disk tends to converge as a result of its viscous evolu¬ 
tion. Setting the left-hand side of equation (8) to zero 
one immediately obtains a steady-state solution for Fj 
in a simple form (Rafikov 2013c) 

Fj{l) = Fj^o + Fj,il = Fj,o + Ml, (10) 

where M = dFj/dl = Fj i is constant. This is different 
from a standard constant M accretion disk solution for 
which M = Stti/S and it follows from equation (4) that 

Fj = Ml = MVLr'^. (11) 

This expression reduces to solution (10) only when Tyo = 
0 and the angular momentum flux at the disk center is 
zero (i.e. no injection of the angular momentum at the 
center by e.g. binary torques). In a cireumbinary disk 
Fj,o 7^ 0 because binary torques inject non-zero angular 
momentum at the center, as / —7 0. As a result, disk 
structure differs from that of a standard constant M disk, 
even though in both cases M is radially constant. This is 
an important distinction that is not always appreciated. 

Even when M = 0 cireumbinary disk has non-zero E 
and dissipation rate, simply because Fj^q ^ 0. In this 
limit of non-accreting binary (Mh ss 0) Fj is radially 
constant, which is illustrated in Figure 1. This is a sit¬ 
uation that we will mainly encounter in this work, see 
§ 6 . 2 . 

All these statements also apply to the quasi-stationary 
disk evolution, in regions where the viscous time is 
shorter than the global evolution timescale of the disk. 
In particular, if the disk is supplied with mass from the 
outside then a solution (10) will develop in its mner parts, 
because the viscous evolution time is short there. This 
solution will have its Eyo and M slowly changing in time 
on the global disk evolution timescale However, Fj 
will still exhibit roughly linear dependence on I out to 
the radius, at which viscous timescale ti, is of order the 
disk lifetime (Rafikov 2013c). 

A remarkable property of the steady-state solution (10) 
is that the behavior of Fj is completely independent of 
the details of physical processes operating in the disk 
— origin of viscous stresses, thermal state of the disk, 
and so on (Rafikov 2013c). Given this solution one can 
trivially obtain the radial scalings of the disk properties 
by solving simple algebraic equations as described in the 
next section. This justifies the introduction of Fj-l vari¬ 
ables and makes interpretation of our results particularly 
transparent. Previous treatments of the cireumbinary 
disk evolution in terms of E and r (Ivanov et al. 1999) 
offered less straightforward interpretation. 

4. CIRCUMBINARY DISKS AROUND STELLAR BINARIES. 

Discussion in the previous section makes it straightfor¬ 
ward that the structure of the cireumbinary disk is best 
parametrized via the angular momentum flux Fj. This 
is in contrast to the conventional accretion disks with¬ 
out angular momentum sources, the properties of which 
are specified via the mass accretion rate M. In the cir- 
cumbinary case such parametrization would be ambigu¬ 
ous, which is easy to understand by looking at the steady 


state solution (10): Fj can vary broadly for a given M 
depending on the value of Fj o- 

To compute detailed properties of cireumbinary disks 
for a given angular momentum flux Fj at some specific 
radius r one needs to specify the thermal structure of the 
disk, since the viscosity depends on the disk temperature, 
see equation (2). Following standard approach, in §4.1 
we calculate midplane disk temperature T by considering 
vertical energy transport in the disk and relating T to E 
and one other disk characteristic. This is a well-known 
exercise in the case of a standard, constant M disk, in 
which T is related to E and M, allowing one to infer 
radial profiles of both E and T. 

By using Fj rather than M as a global disk characteris¬ 
tic, Rafikov (2013c) computed radial behavior of various 
disk properties as a function of Fj in hot and luminous 
cireumbinary disks around supermassive black hole bina¬ 
ries. In this work we extend this approach to circumbi- 
nary disks around stellar binaries. 


4.1. Disk thermodynamics: heating sourees 

Cireumbinary disk around a stellar binary receives en¬ 
ergy from three main sources. First, there is viscous 
dissipation within the disk at the rate dEy/dr per unit 
time and radius given by 


dEy dfl 3 Ejfl 

dr ^ dr 2 r 


( 12 ) 


The associated one-sided energy flux escaping from a disk 
surface is 


1 dEy 3 FM , , 

and is independent of the details of the disk structure 
as long as the value of Fj is specified. In the limit of 
a standard constant M disk with Fj given by equation 
(II) one finds the familiar expression Fy = (3/87r)MD^ 
(Shakura & Sunyaev 1973). 

Second, the disk is also illuminated by the central bi¬ 
nary, which irradiates its surface at a grazing incidence 
angle C ^ 1- Chiang & Goldreich (1997) have demon¬ 
strated that centrally irradiated disks develop a distinct 
two-layer structure: an outer “superheated layer” inter¬ 
cepts direct stellar radiation and re-radiates one half of 
it towards the inner, midplane region, which comprises 
most of the disk mass. The irradiation energy flux re¬ 
ceived by a disk surface is thus^ 




1 Ly 

2 dTrr^ 


c, 


(14) 


where Ly is the combined luminosity of the central stars. 
This expression ignores the time-dependence of irradia¬ 
tion caused by the orbital motion of the central binary 
(Clanton 2013). 

If stellar radiation is intercepted far from the binary 
at a height r]h(r) above the disk midplane then the inci¬ 
dence angle ( can be approximated as 


c 


dr \ r J 


(15) 


This expression assumes that central cavity in the disk allows 
an unobscured view of full stellar disks of both stars from any point 
on the disk surface. 
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Typically ry Ri 3 — 4, and varies with r very weakly (Chi- 
ang & Goldreich 1997). 

There is also a third energy source — shock damping of 
the density waves excited by the binary potential, which 
effect the angular momentum exchange between the bi¬ 
nary and the disk. The one-sided contribution of this 
tidal heating per unit area is 

-^tid = ^ AS, (16) 

where is the angular frequency of the binary and the 
torque A due the binary is given by equation (3). 

4 . 2 . Disk thermodynamics: radiation transport 

We characterize the thermal state of the disk by its 
midplane temperature T and take it to represent the 
characteristic temperature across the whole vertical span 
of the disk, which is of course just an approximation. 
To compute T given the energy fluxes iFv, Jirr and J^tid 
we assume that energy is transported in the vertical di¬ 
rection by radiation and characterize transport by the 
midplane-to-surface optical depth r = (1/2 )Ek. Here k 
is the opacity, which in general is a function of T and 
(midplane) gas density p, which can be easily related to 
S or Fj. 

In the optically thin limit r 1 Kirchhoff’s law and 
energy balance suggest that 

K. Fy F Ftid + rFi-,^. (17) 

Note that technically the values of r in the left and right 
hand sides of this relation are different: the former is 
for the disk opacity at temperature T, while the latter 
characterizes optical depth of the disk to the radiation of 
super-heated dust grains at the disk surface (Chiang & 
Goldreich 1997). The two are somewhat different (disk 
surface is hotter) but for simplicity we will neglect this 
distinction when calculating the midplane temperature. 

In the optically thick case r ^ 1 standard solution for 
radiation transfer in a plane-parallel atmosphere (Calvet 
et al. 1991) yields 

(TT^ Ri —T [FV F F tid) -t- .Tdrr- ( 18 ) 

This expression assumes that viscous energy release oc¬ 
curs close to midplane. Note that in equations (17) 
and (18) we do not differentiate between the Planck and 
Rosseland mean opacities. 

To cover both opacity limits (and interpolate in the t 
1 regime) we use the following expression for T, which 
provides smooth transition between the two regimes: 


aT‘^ = 

“ /(t") .+-741(1)+ -771.., 

(19) 

/(i-)- 

3 1 

(20) 


It is clear that in the limits r <C 1 and t ^ 1 this expres¬ 
sion reduces to equations (17) and (18) respectively. If 
vertical energy transport is affected not by radiation but 
by convection then the relation (19) remains the same 
but the form of f(r) changes in the optically thick case 
(Rahkov 2007). 

Equations (2), (4), (19) form a closed set of equa¬ 
tions uniquely determining E and T for given Fj and r. 


This fully specifies disk properties, including circumbi- 
nary disks which are not in steady state. 


4 . 3 . Irradiation-dominated disk regions. 

We will now describe the disk structure in a partic¬ 
ular important limit when the midplane temperature 
is determined primarily by stellar irradiation, i.e. when 
•^irr 3> /(r) {Fy -I- J^tid), See equation (19). As we will 
see later, this limit is naturally realized in the outer, 
cold parts of the circumbinary disk. Galculation of the 
temperature profile in this regime is essentially identi¬ 
cal to that in the single star case, but we show it for 
completeness as the result is then used in §5. 

In the irradiation-dominated regime the midplane tem¬ 
perature is 



[C(7-)Tc1 

V - ) 

87rCT 


-I 1/4 


.-1/2 


( 21 ) 


Note that the incidence angle ( is itself a function of T, 
according to equation (15). Plugging the expression (15) 
for ^ into equation (21) and solving it, we find (Chiang 
& Goldreich 1997; Rafikov 2006) 


T(r) = 


1? Ac Y 
7 dTTCT / 


;120 K 


viK 

P2Mc,1 


(GMy) 
-\ 1/7 


1/7 


.-3/7 


( 22 ) 


,.-3/7 


In the numerical estimate we use a shorthand notation 
ri = r/AU, r ]3 = p/S, Lcp = L^/Lq, 
and /i is normalized by the H 2 molecular weight. The 
temperature distribution in (22) is independent of Fj. 

The aspect ratio of an externally irradiated disk is 
given by 


h 

r 


a 

fk \ 

7 



At:u(GMc)‘^ 


1/7 


,2/7 


0.024 


V^Lc,! 


1/7 


2/7 


(23) 


SO that the angle at which stellar radiation impinges on 
the disk surface is 


C^(2/7)(Vr). 

From equations (2) and (4) we also find 




■n Lc 

7 Aira 


1 n3 -2 Pj,38 

' 10 g cm - 

a_2 


(GM,) 

(kipf 

viLl, 


1/7 


.-11/7 


1/7 


(24) 


(25) 


-11/7 


where q :_2 = a/lO”^ and Fyag = Fj/(10^®erg). The 
characteristic value of Fj used here is motivated in §5.2, 
see equation (30). 

Note that within the framework of our approximations 
(i.e. not differentiating between the opacity for the ra¬ 
diation of the disk midplane and the superheated outer 
























6 


layer) the properties of passive regions of the disk are 
independent of the opacity behavior and optical depth. 


5. CIRCUMBINARY DISK EVOLUTION: ANALYTICAL 
PICTURE 


Previously obtained results allow us to characterize 
structure and evolution of a circumbinary disk. We as¬ 
sume that most of the disk mass is initially deposited at 
large radii (tens of AU) by the collapse of a centrifugally 
supported envelope. We will also focus on a particular 
case of a disk with no or weak mass inflow into the cavity, 
so that Mb Ri 0. 

After initial mass deposition at some radius vq the disk 
will evolve under the action of viscous stresses. Some of 
the mass will spread inwards towards the binary until it 
reaches semi-major axis ^ (2 — 3)ab. At this point the 
binary torque will become strong enough at imparting 
angular momentum to the disk fluid that the gas inflow 
will be stopped (or at least strongly suppressed) and a 
cavity will form at the disk center. At the same time 
most of the mass remains in the outer disk (see §5.2), 
which spreads out viscously. As a result, the outer ra¬ 
dius of the disk will grow in time. Provided that Mb ~ 0 
the mass of the disk will be conserved during this evo¬ 
lution while its angular momentum will increase because 
of the binary torque. This distinguishes circumbinary 
disks form the conventional circumstellar disks (Shakura 
& Sunyaev 1973), which evolve preserving their total an¬ 
gular momentum (in the absence of external torques) but 
losing mass to accretion onto the central object. 

Characteristic time on which viscous evolution oc¬ 
curs and Fj, E and other disk properties change at 
some radius r is the viscous time t^{r) = r^jv = 
a~^i}r‘^{fi/kBT{r)). Global disk evolution is set by at 
the radius where most of the mass is concentrated, i.e. 
in the outer disk. This is the region where the midplane 
temperature is determined primarily by central irradia¬ 
tion, so that we can use the results of §4.3 to characterize 
disk temperature behavior. Using equations (2) and (22) 
we find in this case 


tu{r)=i 
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This formula demonstrates that ti, goes down as r de¬ 
creases, implying that the inner regions of the disk will 
tend to evolve faster and attain a quasi-steady state de¬ 
scribed in 3.2. Then, because of our assumption Mb = 0 
and according to the reasoning presented §3.2 the disk 
will tend to converge to Fj{r) = const state, as confirmed 
by our numerical calculations in §6.2. 
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This relation is expected to be accurate only for t ^ 
t^{rQ). Note that it predicts r-mfi to be independent of 
the disk mass (or Fj), which is a consequence of our 
natural assumption that irradiation fully determines disk 
temperature in the outer regions. 


5.2. Characteristic value of Fj. 

We now motivate the characteristic value of Fj = 10^® 
ergs that has been adopted in our equation (25) for an 
externally irradiated disk around a stellar mass binary. 
We do this by relating Fj to the disk mass Md{r) en¬ 
closed within some radius r. Using definition (4) and 
the conventional a-prescription for viscosity i/, one finds 
E = {3Trar^Cg)~^Fj. Integrating this over the radial 
span of the disk with 27rr weighting one finds 


Md{r) 


f dr' 
3a J Cg(r') r' 


(29) 


Note that we dropped the lower integration limit in this 
expression, which is ta. This is justified for any disk 
with r TA in which the outer regions contain most of 
the mass. As equation (29) shows, in a constant Fj disk 
this is a direct consequence of T falling with radius. 

Equation (29) allows us to relate Fj in a constant- 
Fj disk to the total disk mass Md and its outer radius 
r. Beyond 10 AU we expect Cg to be determined by 
external irradiation, so substituting T{r) from (22) into 
(29) and performing the integral we express Fj as 


Fj = - 


9 Md 

(■n 


14^3/7^ 

' d 

\ 7 4770- y 

(GMc) 

10®®erg Md,- 2 a -2 

'vlLh ( 

l4Mc,i V 


1/7 


(30) 


n 1/7 


where Md -2 = Md/(10 ^M©). This estimate justifies 
the characteristic value of Fj = 10®® ergs used in equa¬ 
tion (25) for a disk with a typical mass 10 “^Mq and size 
r = 50 AU. 


5.1. Disk expansion 

We now take a closer look at how the disk expands 
in the long run, after the age of the system t exceeds 
the viscous time t^{rQ) at the initial mass deposition ra¬ 
dius tq. We characterize disk expansion by the radius 
of influence rinfl(t), which is the radius at which viscous 
timescale equals the age of the system t. At late times 
Tinfl can be obtained by inverting the relation (26) so 


5.3. Evolution of Fj. 

We can now predict the late time (for t > t^{r[j)) be¬ 
havior of Fj in a constant Fj disk. As Md remains con¬ 
stant because of our assumption Mb Ri O, while the disk 
size increases due to viscous spreading, equation (30) im¬ 
plies that Fj will go down with time. Approximating 
the disk size at a given moment of time t with rinfl(t), 
we substitute the expression (27) for r in equation (30) 
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to find 

(31) 


(32) 

This formula determines the scalings of the characteristic 
Fj in a constant-Fj disk with the disk mass and the 
age of the system t, which will be subsequently verified in 
§6.2. It can be used to determine the evolution of E(r, t) 
in the outer disk dominated by central irradiation with 
the aid of equation (25). 

Note that equation (31) does not contain any details 
of the opacity behavior in the outer disk. This is a direct 
consequence of our simplified treatment of the radiation 
transfer, in which we took the dust opacity to be the 
roughly the same for both the midplane temperature T 
and the temperature of the superheated surface layers of 
the disk (recognizing this distinction would would give 
rise to only small difference with our results). 

6. CIRCUMBINARY DISK EVOLUTION: NUMERICAL 
RESULTS. 

Now we describe the results of our ID numerical calcu¬ 
lations of the circumbinary disk evolution based on equa¬ 
tion (1). In doing this we provide close comparison with 
the analytical predictions outlined in §5. To better illus¬ 
trate the differences between the circumbinary and stan¬ 
dard constant M circumstellar disks, we also performed 
evolutionary calculations for disks orbiting a single star 
with the same total mass Me and no torque injection 
at the center. Our numerical approach and parameters 
of our simulations are described in Appendix A, which 
should be referred for details. 

6.1. Viscous expansion of the disk 

First, we explore the viscous expansion of the disk seen 
in our simulations. As a proxy for the outer radius of 
the disk we use the radius tm at which the “effective 
disk mass” E(r)r^ attains its maximum value. Figure 2 
shows a plot of rM{t) obtained from our simulations for 
different disk masses. 

At early times, for t < t^fro), one naturally has tm ~ 
vq reflecting the radius of initial mass deposition. Later 
on, for t > ti,{ro), the disk viscously expands and rM(t) 
grows. By fitting our numerical results we find that at 
late times 
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where rinfl(t) is given by equation (27). 

Note that late-time tm {t) is only weakly dependent of 
Md in our simulations. This is in agreement with our an¬ 
alytical expectations (see §5.1), because disk expansion 
is set by the physics of the outer, externally-irradiated 
disk regions, in which both T and v are independent of 
Md, see equation (22). 

6.2. Evolution of Fj distribution. 



Fig. 2.— Semianalytic fit to the evolution of 7’infi(f). Solid lines 
are the results of the numerical calculations, and the dashed line is 
tmIF) given by equation (33) and based on analytic theory (equa¬ 
tion (27)), which fits numerical results at late times. Different 
colors correspond to different disk masses: black Md = O.lMc, red 
Md = O.OSMc, and green Md = O.OlMc. 

We now explore the evolution of Fj distribution in our 
numerical models. Figure 3 shows Fj(r,t) at different 
moments of time for a disk with Md = O.lMc- 

Initially the disk evolves on a timescale considerably 
shorter (by an order of magnitude) than the viscous time 
tv = T%lv at To, which is in the irradiated part of the disk. 
Figure 3 shows that already at 5 x 10® yr « 0.01tj/(ro) 
(black curve) viscous mass inflow reaches the central bi¬ 
nary, where it is stopped by the tidal torque, giving rise 
to a plateau in the radial distribution of Fj at r of or¬ 
der several AU. Such rapid evolution is an artefact of the 
strong radial gradients of E associated with its initial 
narrowly-peaked distribution in the form of a ring, see 
equation (Al). Sharp drop of Fj at small radii is caused 
by the disk truncation by the binary torques and forma¬ 
tion of the central cavity. At larger separations r ^ ro 
relaxation to Fj = const state has not yet been achieved, 
and mass flows both inward and outward from the r ^ ro 
region. 

At t = 7 X 10"^ yr r; 0.16F(ro) (blue curve) radial dis¬ 
tribution of Fj becomes consistent with being constant 
over more than an order of magnitude in radius. This 
implies that M in this part of the disk is very small. 
At later times viscous expansion of the disk extends the 
outer disk edge by about an order of magnitude, gener¬ 
ally preserving radially constant Fj over the large span 
of the disk. 

At the same time, the amplitude of Fj in the bulk of 
the disk (the height of the plateau) constantly decreases 
at late time, in agreement with theoretical expectations 
(§5.3). To quantify this behavior we measure Fj* — the 
height of the plateau of the Fj(r,t) distribution in our 
simulations — for three different disk masses at different 
moments of time and plot it in Figure 5. It is clear that at 
late stages, t >tv(ro), Fj decays roughly as a power law 
in time. Also, more massive disks feature higher values 
of Fj. 

By fitting these numerical results we found that at late 
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Fig. 3.— Evolution of Fj{r, t) in a circumbinary disk with = 
O.lMc. Different curves correspond to different epochs: black — 
5.0X10^ yrs, blue — 5.0x10^ yrs, green — 5.0x10^ yrs, red — 5.0x 
10® yrs. At time t = 0 mass is distributed in a narrow Gaussian 
ring at radius ro = 20 AU, see equation (Al). It subsequently 
spreads viscously both inward and outward of rp. Note the overall 
expansion of the disk with time and formation of the central cavity. 



t (yrs) 

Fig. 4.— Time evolution of Fj^{t) — value of Fj at the plateau 
portion of the Fj(r) distribution in our simulations of circumbi¬ 
nary disks around non-accreting central binary. Different colors 
correspond to different disk masses: = O.lMc (black), 0.05Mc 

(red), and O.OlMc (green). Solid curves are the results of our nu¬ 
merical calculations, while dashed lines represent our analytical fit 
given by equations (31) & (34). 

times 

Ff{t) « 0.45Fj(t), (34) 

where Fj{t) is given by equation (31). This verifies and 
calibrates our simple scaling (31) for the behavior of Fj 
presented in §5.3, and motivates the use of our analytical 
results in other circumbinary disk applications. 

Evolution of Fj distribution shown in Figure 3 should 
be contrasted with Fj{r,t) behavior in the case of a cir- 
cumstellar disk around a single star without angular mo¬ 


r/ro 

0.01 0.1 1 10 100 



r (AU) 

Fig. 5.— Same as Figure 3 but for a circumstellar disk of 
Md = O.lMc around a single star with the same total central mass 
as in circumbinary case. Different colors correspond to the same 
moments of time as in Figure 3: black — 5.0 x 10^ yrs, blue — 
5.0 X 10"^ yrs, green — 5.0 x 10® yrs, red — 5.0 x 10° yrs. How¬ 
ever, the structure of the radial distribution of Fj is very different 
compared to the case of a circumbinary disk. 

mentum injection at the center. Figure 5 displays how Fj 
evolves in the latter case, keeping the same central mass 
Me and all disk parameters. Profiles of Fj are plotted 
at the same moments of time as in Figure 3. One can 
see a dramatic difference in Fj behavior compared to the 
circumbinary case: in circumstellar disk Fj(r) never be¬ 
comes flat. Instead, Fj distribution rapidly converges to 
a standard constant M solution (11) with Fj,q = 0, i.e. 
Fj(r) oc 

The amplitude of Fj in Figure 5 rapidly drops in time 
for two reasons. First, analogous to the circumbinary 
case, the disk expands, lowering E and, consequently, re¬ 
ducing Fj (this effect operates in circumbinary disks as 
well). Second, ongoing mass accretion onto the central 
object (absent in the circumbinary disk with = 0) 
reduces disk mass on global viscous timescale and addi¬ 
tionally lowers Fj. Lower value of Fj in a circumstellar 
disk, especially in its inner regions, diminishes the role 
of viscous heating in determining its thermal state com¬ 
pared to the case of a circumbinary disk. 

6.3. Evolution of disk properties. 

We will now look at how other basic disk properties 
evolve in our numerical models. Figures 6-8 summarize 
our results for three different disk masses, Mj = O.lMc, 
0.05Mc, and O.OlMc, correspondingly. Each figure shows 
the surface density S(r), midplane temperature T(r), 
and optical depth r(r) at the same moments of time as 
in Figure 3. 

Examination of these figures reveals several features 
common to circumbinary disks of different masses. The 
overall shapes of S(r) curves are similar in different plots, 
despite the overall scales being different at the same mo¬ 
ments of time: higher E for more massive disks, and vice 
versa. Surface density increases towards smaller radii 
and exhibits a clear signature of a central cavity that is 
placed at 2ab by design (by choosing a proper value of 
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Fig. 6.— Evolution of the surface density (top), midplane tem¬ 
perature (middle), and optical depth (bottom) in our numerical 
models of a high-mass circumbinary disk with = O.lMc with¬ 
out mass accretion at the center, Mj, ^ 0. System with the follow¬ 
ing parameters is shown: equal mass (q = 1.0) central binary of 
total mass Me = Mq and luminosity Lc = 2 Lq,. For this disk mass 
we set / = 2 X 10~^ in the torque density prescription (3) for the 
central cavity to have a radius Vq ~ 2ab, in agreement with hydro 
simulations. Different colors correspond to different times: black 
5.0 X 10^ yrs, blue 5.0 x 10^ yrs, green 5.0 x 10^ yrs, red 5.0 x 10® 
yrs. The dotted lines in two upper panels correspond to the ana- 
l 3 d.ic scalings for E(r) and T(r) for the irradiation-dominated case, 
equations (22) and (25). Note the relatively slow evolution of S(r) 
caused by the assumed lack of accretion onto the binary. 

the coefficient / in our torque density prescription (3)), 
indicated in each Figure. Peak values of S('r) at each 
moment of time scale roughly linearly with used in 
the calculation. 

At the same time, temperature profiles for different 
disk masses vary less significantly, with the difference 
in peak values of ^ 2 between the high- and low-mass 
disk models (M^ = O.lMc and O.OlMc, respectively). 
Moreover, beyond 10 — 30 AU the midplane temperature 
converges to the prediction (22) at all times, because 
central irradiation becomes the dominant factor setting 
T (see §8.1). For the same reason E(r) on these scales 
obeys equation (25) at late times (dotted lines), see the 
red curve at 5 Myrs in top panels of Figures 6-8. 

Viscous dissipation plays important role in setting mid¬ 
plane temperature of the inner disk regions. It domi¬ 
nates over irradiation for r < 10 — 30 AU for Md = 
(0.05 — 0.1)Mc, while for Md = O.OlMc its effect is con¬ 
fined to within several AU from the binary. High-mass 
disks can be warmed up by viscous and tidal dissipation 
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Fig. 7.— Same as Figure 6 but for Md = O.OSMc and / = 
1.5 X 10“^. 

to 2 X 10^ K in their inner regions, while in low mass 
disks T does not exceed 10^ K. 

Optical depth of the disk (lower panels of Figures 6-8) 
exhibits complicated structure related to sublimation of 
different grain species responsible for the opacity of the 
disk (Zhu et al. 2009). These opacity transitions are one 
of the reasons for the fine structure visible in S(r) and 
T(r) profiles. They become less numerous at lower disk 
masses since viscous heating and tidal dissipation do not 
raise midplane temperature to very high values necessary 
for sublimating more refractory grains in this case. 

All our numerical models become optically thin beyond 
40-100 AU, depending on Md- There T{r) may differ 
slightly from the prediction (22), see the discussion after 
equation (17). However, the correction — a weak power 
of the ratio of grain emissivity at the peak wavelengths of 
the blackbody emission of the superheated layer and the 
disk midplane (Chiang & Goldreich 1997), set to be the 
same in our simple approximation — is not very different 
from unity and varies with r very weakly to change our 
results significantly. 

Now we provide comparison with the case of a circum- 
stellar disk without external sources of angular momen¬ 
tum. In Figure 9 we show the evolution of Md = 0.05Mc 
circumstellar disk; it should be compared with Figure 
7 which illustrates the evolution of a circumbinary disk 
with Mf, ss 0 of the same mass. In computing the evolu¬ 
tion of a circumstellar disk we simply set A = 0, keeping 
everything else the same as in the calculation of Figure 
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Fig. 8.— Same as Figure 6 but for = O.OlMc and / = 
1.0 X 10-3. 

7. 

Lack of central cavity is not the only obvious feature 
of the circumstellar disk. Comparing Figures 9 and 7 we 
see that E drops much faster with time in a circumstel¬ 
lar disk. Also, the slope of ^{r) profile is shallower in 
the circumstellar case, see the red curve in top panel. It 
obeys E oc which can be obtained from equation 

(25) by substituting Fj oc as appropriate for a cir¬ 
cumstellar disk. These differences are easily understood 
in terms of the very different evolution of the Fj{r, t) dis¬ 
tributions shown in Figures 3 and 5: in the circumstellar 
case Fj (and, subsequently, E) is much lower in the inner 
disk and drops with time more rapidly due to ongoing ac¬ 
cretion onto the central object. Because of that at late 
time (t ^ 5 Myr) T(r) in circumstellar disk is set entirely 
by irradiation outside of ^ 0.5 AU (compared to ~ 10 
AU in a circumbinary case). 

Based on that we conclude that circumbinary disks are 
very different from their circumstellar analogues: they 
tend to remain more massive and maintain high levels of 
viscous dissipation, which keeps them hotter for longer. 

7. SPECTRA OF CIRCUMBINARY DISKS. 

Given significant differences in properties of circumbi¬ 
nary and circumstellar protoplanetary disks outlined in 
previous section one should expect their spectral energy 
distributions (SEDs) to also differ. We describe the de¬ 
tails of our SED calculations in Appendix B and provide 
comparison between SEDs of two types of disks in this 



1 10 100 
F (AU) 

Fig. 9.— Same as Figure 7 but for a circumstellar disk with non¬ 
zero accretion onto the star but no external angular momentum 
injection. Dotted line in the upper panel illustrates E(r) oc j-—15/14 
behavior typical for irradiated disks. See §6.3 for details. 

section. 

Figure (10) shows the SEDs of the circumbinary (top 
panel) and circumstellar (middle panel) disks with start¬ 
ing mass Md = 0.05A/c at t = 5 Myr. In these plots 
we separately show the contributions from the stellar 
blackbody (dotted), superheated dust layer (long-short 
dashed), and the bulk of the disk (dashed). The last 
component is also split into tidal (only for circumbinary 
case), viscous, and irradiation contributions. Each of 
these is evaluated by setting the other heating sources to 
zero®. 

The bottom panel of Eigure (10) provides a direct com¬ 
parison of the total SEDs obtained for the two disks. Its 
inspection highlights three important differences. First, 
the circumstellar disk is a stronger emitter at A ^ 3/rm 
simply because it gets closer to the star and reaches 
higher temperatures. Presence of the central cavity sup¬ 
presses the near-IR emission of the circumbinary disk. 

Second, circumbinary SED exhibits a bump at A ~ 
lO^m, which is absent in circumstellar SED. This fea¬ 
ture arises because of the tidal dissipation in a circumbi¬ 
nary disk, which dominates heating of its inner parts (see 
§8.1). Note that in our SED calculation the significance 
of the bump is reduced by the high level of the super- 

® I.e. tidal contribution is obtained by setting y-j) = = 

0 in equations (Bl) and (B3); other contributions are obtained 
similarly. 
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Fig. 10. — SEDs of a circumbinary (top panel) and circumstellar 
(middle panel) disk with = 0.05Mc. Different colors corre¬ 
spond to different contributions to the SED: solid black (total), 
dotted black (stellar), violet (superheated grains), and red (disk). 
The disk contribution is further divided into its components: blue 
(tidal heating), cyan (irradiation), and green (viscous heating). 
The middle panel, corresponding to the SED of a circumstellar 
disk, has no tidal component. The bottom panel compares the to¬ 
tal circumbinary SED (solid black line) with the total circumstellar 
SED (dotted black line). 

heated dust emission, which dominates because we ne¬ 
glect the fact that upper layers of the disk should become 
optically thin to stellar radiation at large separations. If 
we accounted for this effect, the superheated dust contri¬ 
bution would be lower, giving rise to larger difference be¬ 
tween the circumbinary and circumstellar SEDs around 
lO^m. 

Third, circumbinary disk provides more emission in the 
Rayleigh-Jeans tail. This is because at t = 5 Myr this 
disk has more mass than its circumstellar counterpart, 
which has already lost « 82% of its original mass to 
accretion onto the central star by that time. 

The distinctive bump in the SED at A ~ lO/rm may 
serve as a circumbinary disk signature that could be used 
to identify such systems in observations, when the bina- 
rity of the central source is difficult to determine. Such 
SED-based technique using features in quasar spectra 
has been previously suggested as a way of finding binary 
supermassive black holes in centers of galaxies (Rafikov 
2013c; Yan et al. 2015). In practice, in protostellar case 
such identihcation may be complicated by the presence 
of a strong silicate dust resonance in the disk emission at 
9.7/im (Mathis 1990; Rafikov & De Colle 2006) (we do 
not account for it here), which could lead to ambiguity 
in interpretation. 



Fig. 11.— Aspect ratio of O.OSMq circumbinary disk. Different 
colors correspond to different times: black 5.0 X 10^ yrs, blue 5.0 X 
10^ yrs, green 5.0 X 10® yrs, red 5.0 X 10® yrs. The disk is self- 
shadowed at 1-20 AU. 

7.1. Disk shadowing 

Our detailed calculations of the thermal structure of 
circumbinary disks in §6.3 reveal an interesting effect, 
which is illustrated in Figure 11. In this Figure we plot 
the disk aspect ratio h/r as a function of radius at dif¬ 
ferent moments of time. The outer parts of the disk are 
dominated by central irradiation and h/r behavior there 
is well described by equation (23), resulting in a flared 
disk structure. However, at r < 10 AU situation changes 
and h/r stays roughly constant or even increases with 
decreasing r. In our calculations this pattern persists at 
all times (except the earliest epochs affected by the ini¬ 
tial conditions) and implies that between 1 and 10 AU 
circumbinary disk should have non-flared structure. 

This behavior is caused by the vigorous tidal and vis¬ 
cous heating in the inner disk, as we demonstrate in de¬ 
tail in §8.1. This energy release dramatically increases 
midplane temperature of the disk and boosts h/r. This 
feature is unique for circumbinary disks (i.e. it is not 
present in the circumstellar disks) because of the tidal 
heating and peculiar Fj{r) behavior endemic for such 
disks. 

Deviation oi h/r behavior from equation (23) means 
that our procedure for calculating (which assumes C, 
to be given by equation (24), see §A) becomes inaccurate 
at r < 10 AU. Moreover, in regions where h/r is lower 
than in some interior part of the disk, that inner por¬ 
tion of the disk blocks stellar irradiation from reaching 
the disk surface and J-i„ —> 0, resulting in self-shadowing 
of the disk. Accounting for this effect is greatly compli¬ 
cated by its global nature of the effect. Accurate descrip¬ 
tion of the self-shadowing may necessitate direct radia¬ 
tion transfer calculations. 

Both our evolutionary and SED calculations presently 
include reprocessed stellar radiation emitted from disk 
regions that are in fact shadowed. This is not a serious 
issue in the inner disk (< 10 AU) simply because there 
-Firr is subdominant relative to both Ty and J^tid- Thus, 
retaining irradiation contribution in the energy budget 
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produces negligible effect in this region. 

More worrisome could be the fact that the inner disk 
regions may cast shadow over significant portion of the 
disk at intermediate separations where we currently con¬ 
sider central irradiation to dominate over the viscous and 
tidal heating. For example, focusing on 0.5 Myr curve 
in Figure 11 we see that aspect ratio peaks at the level 
of h/r Ri 0.07 at r « 2 AU. Even though Jirr starts to 
dominate h/r behavior around 10 AU, aspect ratio again 
reaches 0.07 in the irradiation-dominated outer disk only 
around 20 — 30 AU. This implies that the geometrically- 
thick parts of the disk at 2 AU cast shadow out to 
20 — 30 AU at this moment of time. This is likely to 
lower temperature in this part of the disk and reduce 
irradiation contribution® to the SED coming from that 
region. While we neglect self-shadowing in this work, 
future studies should account for its effect on the disk 
structure and observational appearance. 

8. DISCUSSION. 

We now discuss some additional aspect of the circumbi- 
nary disk evolution which follow from the results pre¬ 
sented in §6. 

8.1. Role of different heating terms. 

Here we explore the role of different heating sources 
in determining the thermal state of a circumbinary disk 
at different radii. Figure 12 shows the runs of viscous 
J7u(r), tidal -Ftid(?'), and irradiation Jdrr('c) heating at 
different moments of time for Md = O.OSMc disk. Several 
important conclusions can be drawn from this plot. 

First, on scales r < 2 AU dissipation of the density 
waves launched by the binary provides the most impor¬ 
tant heating source. At its peak (around 0.6 AU) J^tid 
exceeds Ty by about an order of magnitude at all times. 
To understand this behavior we will use the fact that 
viscous time in the inner disk is much shorter than the 
global evolution timescale, so that a quasi-steady state 
must develop there. In this limit the binary torque near 
the disk edge is balanced by the gradient of the viscous 
angular momentum flux, so that dFj/dr = 27 rrEA. We 
can then use this relation and equations (13) and (16) to 
write 

J-tid 2 ainfjUb-H ^ 2 91nfj / r 
Fy 3 cllnr 11 3 cllnr 

where the last approximation follows from the fact that 
^ everywhere in the disk. 

Looking at the Figure 3 we see that dlnFj/d\nr is 
very large near the inner edge of the disk. This, coupled 
with the fact that (r/a^)®/^ ^ 1, explains the major role 
of Ffid in heating the disk within ~ 2 AU. Outside this 
region dlnFj/dhir rapidly decreases {Fj distribution 
flattens, see Figure 3), and viscous dissipation starts to 
dominate over the tidal heating. 

The dominant role of the tidal heating in the inner part 
of the disk makes it important to understand the detailed 
radial structure of the density wave dissipation in the 
disk, since J^tid is directly connected to A(r). Our simple 
prescription (3) ignores both the recent developments in 

® This would additionally accentuate lO/xm bump in the cir¬ 
cumbinary SED discussed in §7. 



Fig. 12.— Relative contributions of different heating terms for a 
0.05Mq circumbinary disk: viscous heating (dashed), tidal dissi¬ 
pation (dotted), and irradiation (solid black). From top to bottom 
the times are: 5.0 X 10® yr (black), 5.0 X 10^ yr (blue), 5.0 X 10® yr 
(green), 5.0 X 10® yr (red). Irradiation does not evolve with time 
and, at late times, viscous heating is dwarfed by irradiation for all 
radii. Tidal heating dominates in the inner disk. 



Fig. 13.— Same as Figure 12 but for a O.O5M0 circumstellar disk 
(note the lack of tidal heating). Viscous heating is less significant 
in the circumstellar case than in the circumbinary case because of 
the pile-up of mass at small radii in the latter case. 

understanding the excitation torque density (Rafikov & 
Petrovich 2012; Petrovich & Rafikov 2012) and the non¬ 
local nature of the density wave damping (Goodman & 
Rafikov 2001; Rafikov 2002), motivating further refine¬ 
ments. 

Second, irradiation by the binary is always the major 
heating source at large separations. This is because Fy oc 
in a constant-Fj disk, while Firr cx; C,{r)/r'^ decays 
much slower with r. Because both viscous heating and 
tidal dissipation scale with Md^ we expect irradiation to 
dominate closer to the star for less massive disks. 

Third, even when > Fy , Ftid tidal and viscous 
dissipation may not be neglected when calculating the 
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midplane disk temperature T. For example, Figure 12 
shows that J^irr ~ J-'u + J^tid around 1.5 AU at t = 5 
Myr. However, Figure 7 demonstrates that T converges 
to the behavior (22) only outside 8 AU. This is because 
J-y and Jyid have a weighting factor r in equation (19) 
determining T, and the disk is optically thick inside 10 
AU. 

In Figure 13 we plot J^rr and in a circumstellar 
disk of the same mass as in Figure 12. Comparing the 
two Figures one can see that in the circumstellar case 
J-y becomes a subdominant heating source faster than in 
the circumbinary case, as a result of accretion onto the 
central star in the former case. Also, the Ty{r) profile 
is shallower in the circumstellar case, because Fj oc 
for a constant-M disk (see §6.2) and equation (13) then 
predicts Fy{r) oc r~^, as opposed to Fy oc in a 

constant-Fj circumbinary disk. 


r < nnfl) is given by 

Fj{r) Ki Fj(ri„fl) + Mb [l{r) - /(ri„fl)]. (36) 

With this relation one can then easily determine the pro¬ 
files of S(r), r(r), etc. using definition (4) in much the 
same way as we derived equation (25). 

This simple recipe provides a prescription for quan¬ 
titative understanding of the effect of non-zero binary 
accretion on evolution of the circumbinary disk proper¬ 
ties. Apparently, since the disk is losing mass, Fj is going 
to be smaller than in the constant-Fj circumbinary disk 
with Mb ~ 0, resulting in lower-E, cooler, and less lu¬ 
minous disk. It is also clear that properties of the disk 
with Mb ^ 0 will always be somewhere in between the 
circumstellar case with no torque (which provides the 
maximum Mb) and the constant-Fj circumbinary disk, 
both of which have been covered in §6. 


8.2. Effect of accretion onto the binary 

In most of this work we assumed that the binary torque 
cuts off disk inflow into the central cavity and Mb « 0. At 
the same time, Artymowicz & Lubow (1996) suggested 
that tidal streams could penetrate the cavity, giving rise 
to some accretion by the binary. Early work (Mac- 
Fadyen & Milosavljevic 2008) suggested that Mb due to 
tidal streams is ^ 10% the accretion rate without cen¬ 
tral torque. However, recent hydrodynamic (Farris et al. 
2014; D’Orazio et al. 2013) and MHD (Shi et al. 2012) 
simulations suggest that Mb may reach (30 — 60)% of the 
corresponding accretion rate without central torque. 

For brevity we do not show numerical disk models with 
Mb 7 ^ 0 (except or the standard constant-M circumstel¬ 
lar case) in this work. However, we outline qualitatively 
how our results should change in this case. 

As long as the outer disk edge remains in the 
irradiation-dominated part of the disk equations (26) and 
(27) remain valid. But equation (29) gets modified for 
two reasons. First, because of accretion at some (gener¬ 
ally time-dependent) rate Mb{t) disk mass is now a de¬ 
creasing function of time and dMd/dt = —Mb- Second, 
Fj{r) is no longer constant with r and instead obeys the 
solution (10), as long as variations of Mb are not faster^ 
than the global disk evolution. For that reason formally 
we can no longer take Fj out of the integral in equation 
(29). 

Nevertheless, it is obvious that this integral is still dom¬ 
inated by the outermost disk regions, where Fj attains 
its maximum value ~ Fj(rinfl). Thus, we can approxi¬ 
mately evaluate the integral at its upper limit and find 
that, up to factors of order unity, equations (30) and 
(31) remain valid as long as we (1) replace in them Fj 
with F,/(rinfl) and (2) consider disk mass to be a func¬ 
tion of time, Md{t), which can be computed once the his¬ 
tory of accretion onto the binary, Mb{i), is known. Once 
the approximate value of Fj(rinfl) is determined and cur¬ 
rent value of Mb is known, the global run of Fj{r) (for 


^ It has to be kept in mind that Mbif is likely to be set not by 
the global disk properties, but by its characteristics near the cavity 
edge, since this is where binary torque is posing a barrier to gas 
inflow. 


8.3. Binary Inspiral due to Disk Coupling 

Throughout this work we assumed binary separation 
to remain fixed at Ob = 0.2 AU. However, in practice, 
binary is constantly losing its angular momentum due to 
tidal coupling with the disk. Here we look at this process 
in some detail. 

The amount of angular momentum that the binary 
possesses is 


Lb 


q 

(1 -b g)2 


{GM^abY^^ 


(37) 


4 X 10®^ g cm^s ^ 


q 

(i + g)2 


M. 


3/2 

C,1 



Binary torque provides a source of the angular momen¬ 
tum for the disk with ” angular momentum power” equal 
to Fj_o ~ Fj(rA), where Fyo is featured in quasi-steady 
solution (10) and ta is the radius, beyond which only 
viscous stresses matter (i.e. binary torque can be ne¬ 
glected). As a result Lb changes at the rate Lb = —Fj,o. 

In our constant-Fj disk with Mb « 0 we have Fyo ~ 
Fj* — the value of Fj at the plateau of its radial dis¬ 
tribution. Using equations (31) and (34) based on our 
numerical results one can easily compute the change of 
Lb accumulated over time t to be 


Z 

AFb(t)«0.45 J Fj{t')dt' PC 0.8Fj{t)t 


9 X 10®* g cm^s * Md,_2aL4^^ 

- T 1 2/13 / , \ 7 '^ 

McpFc,1173 I _t _\ 

Myr; 


M 2 


(38) 


(39) 


In Figure 14 we show the amount of the angular 
momentum lost by the central binary as determined 
from one of our runs with Md = 0.05Mc and stan¬ 
dard binary parameters. Solid curve shows AF = 
2'k J* r'K{r')Y,(r')dr' (angular momentum absorbed by 
the disk, which is equal in amplitude to the angular mo¬ 
mentum lost by the binary) as a function of time. The 
dotted line is the theoretical prediction (39), which fits 
the numerical data reasonably well at late times. The bi¬ 
nary parameters adopted in this calculation (e.g. <? = 1) 
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Fig. 14. — Binary angular momentum loss due to tidal interaction 
with a 0.05Mq circumbinary disk in the case of no central accre¬ 
tion, Mf, 0. Dotted line shows the late-time analytical prediction 
for ALf){t) given by equation (39). 

imply that initially Lf, 10^^ g cm^ s“^. Figure 14 
shows that this amount of angular momentum gets lost 
by the binary due to tidal coupling with the disk already 
by 0.4 Myr. Thus, this particular system would merge 
into a single star within that interval of time. 

Note that our calculation of ALi,{t) is essentially in¬ 
sensitive to the binary parameters — binary simply pro¬ 
vides a barrier for the disk inflow causing mass accu¬ 
mulation at the inner edge. The global evolution of the 
disk is independent of the binary characteristics as long 
as Mb remains the same and Fj in the bulk of the disk 
stays unchanged. For that reason our analytical estimate 
(39) does not involve binary semi-major axis, for exam¬ 
ple. Thus, we could have done the same calculation with 
wider or more massive binary, having higher Lb, and still 
used the curve in Figure 14 to determine whether it will 
merge or not (i.e. whether ALb{t) ever reaches the ini¬ 
tial Lb). Using this Figure we find, in particular, that all 
binaries with initial Lb ^ 6 x 10®^ g cm^s“^ (e.g. more 
massive or more widely separated) would avoid merger 
after being embedded in a circumbinary disk for 5 Myr 
(although their orbits could still shrink appreciably). 

Our estimate (39) has been done for the constant- 
Fj disk arising when Mb = 0, and it is natural to 
ask how it would change if some accretion onto the bi¬ 
nary is allowed. In that case Fj profile would be rep¬ 
resented by equation (36), and the amount of angular 
momentum lost by the binary per unit of time becomes 
Fj{r —>■ 0) « Fj{rinfi) — Mb^nnfl)- Once the evolution 
of Tinfl and Fj(rinfl) is specified as described in §8.2, one 
can again compute the angular momentum lost by the 
binary simply as ALb{t) « f* Fj(r 0, t')dtt'. It is obvi¬ 
ous that a disk with non-zero Mb will absorb less angular 
momentum from the binary than its counterpart around 
a non-accreting binary, simply because Fj{r —>■ 0,t) is 
always higher in the latter case. Thus, accreting binaries 
have a higher chance to survive against orbital inspiral 
caused by the tidal coupling to the circumbinary disk. 



1 10 100 
r (AU) 


Fig. 15. — Location of the dead zone for a circumbinary disk of 
mass 0.05Mq at different moments of time: blue 5.0 X 10^ yrs, 
green 5.0 X 10® yrs, red 5.0 X 10® yrs. Filled diamonds on each 
S(r) and T(r) curve represent the inner and outer edges of the 
dead zone for Scrit = 40gcm“^. Open diamonds illustrate the 
outer boundary of the dead zone for l^crit = lOOgcm”^. Critical 
temperature at which thermal ionization becomes effective is taken 
to be 800 K in this calculation. 

8.4. Dead zone 

Following the work of Martin et al. (2013) we inves¬ 
tigate the possibility of the dead zone (Gammie 1996) 
formation in our circumbinary disk models. We assume 
MRI to be inactive and a dead zone to form provided that 
the midplane temperature T < 800 K and disk surface 
density E > 'Lcrit ~ 100gcm“^ (Gammie 1996; Ume- 
bayashi & Nakano 1981; Zhu et al. 2009). At tempera¬ 
tures this low thermal ionization is inefficient, while the 
surface density constraint guarantees that cosmic ray ion¬ 
ization is not effective either. It has to be kept in mind 
that the exact value of Ecr-it is not well known (Turner 
& Sano 2008; Bai & Stone 2011; Zhu et al. 2009) and 
we consider it as a free parameter in our study. We dis¬ 
play the location and evolution of the dead zone for a 
Mil = O.O5M0 circumbinary disk in Figure 15. At each 
snapshot the inner and outer boundaries of the dead zone 
for Ecrit = 40 g cm“^ are shown as hlled diamonds of 
the corresponding color. For higher Ecrit = 100 g cm“^ 
the outer boundary of the dead zone shifts inward (open 
diamonds). We typically find that the dead zone should 
be present in the circumbinary disk and span several AU 
in radius. This is consistent with the findings of Martin 
et al. (2013) even though they used considerably lower 
values of Ecrit than we do here. 

Because of the uncertainties related to the details of 
the dead zone structure (E^rit, value of the background 
viscosity inside the zone, etc.) our present work neglects 
this possibility and focuses on disk evolution with a uni¬ 
form value of a. But our results shown in Figure 15 do 
call for the development of more accurate models of the 
circumbinary disks that would account for the possible 
existence of the dead zone. 
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Fig. 16.— Evolution of the iceline position rice for circumbinary 
(solid curves) and circumstellar (dashed curves) disks. Different 
colors correspond to different starting disk masses: black = 
O.lMc, red = O.OSMc, and green = O.OlMc. 

8.5. Iceline 

Another important location in the disk is the iceline 
— radius rice at which water vapor (and other volatile 
species) condense into solids. The position of iceline in 
the circumbinary disks has been previously addressed 
in Martin et al. (2013), Clanton (2013), Shadmehri & 
Khajenabi (2015), and we provide comparison with these 
studies in §8.7. Here we mainly focus on the differences 
in the iceline location in the circumbinary and circum¬ 
stellar disks. 

In Figure 16 we have plotted rice as a function of time 
for the three values of M 4 for the circumbinary (solid) 
and circumstellar (dashed) disks. Unlike most of the 
prior works (Hayashi 1981; Martin & Livio 2012), which 
define the iceline simply as a location at which the disk 
temperature reaches ss 170 K, here we use the opacity 
tables from Zhu et al. (2009) to account for the pressure 
dependence of the sublimation temperature. 

Our calculations clearly show that for a given start¬ 
ing Md iceline lies at larger distance in the circumbinary 
disk compared to its circumstellar analogue. For exam¬ 
ple, at 1 Myr the circumbinary disk with Md = O.ObM^ 
has iceline at 6 AU, almost a factor of 3 further than in 
the circumstellar disk of the same mass. This difference 
arises because of the presence of strong tidal heating in 
the inner parts of circumbinary disks, which pushes rice 
out. In addition, viscous heating is also more significant 
in the inner circumbinary disk because of both the differ¬ 
ent Fj{r) distribution and the assumed lack of accretion 
onto the central binary, resulting in higher disk mass at 
late times. 

The decrease of rice with time in the circumbinary case 
is driven by the viscous expansion of the disk, which low¬ 
ers its characteristic Fj. In circumstellar case Fj addi¬ 
tionally decreases due to mass loss to accretion, explain¬ 
ing why the difference in 7’ice between the two cases grows 
with time. This has implications for planet formation, 
which are discussed in §9. 

Also, as expected, more massive disks have icelines lo¬ 


cated further out because of the increased contribution 
of viscous and tidal heating, both of which scale with 
disk surface density or Fj, see equations 13 & 16). 

8.6. Limitations of our model. 

We now comment on the validity of the approximations 
used in this work and outline its limitations (in addition 
to not accounting for the possibility of a dead zone, see 
§8.4). 

Our model ignores the azimuthal structure of the disk 
which should arise due to the rotating potential of the 
binary. However, in the time-averaged sense our axisym- 
metric model should still represent fairly well the actual 
disk structure. This is especially true beyond severalx a;,, 
where the time-varying non-axisymmetric disk structures 
become subdominant. The same is true for the variabil¬ 
ity in time and space of the irradiation flux impinging on 
the disk surface (Clanton 2013; Shadmehri & Khajenabi 
2015; Hodman & Quillen 2015). 

Our treatment of disk thermodynamics has room for 
improvement. Our calculation of the midplane disk tem¬ 
perature ignores the possibility of the disk self-shadowing 
in regions where viscous and tidal heating dominate, 
see §7.1. This motivates future more detailed, higher¬ 
dimensional studies of the radiation transport in the cir¬ 
cumbinary disks that would be able to account for disk 
self-shadowing. 

Our assumption of equality between the temperature 
of the superheated grains and the midplane temperature 
in the irradiated case (see §7) may affect the determi¬ 
nation of T when the disk becomes optically thin. The 
rate at which the disk viscously spreads determines evo¬ 
lution of many of its properties (e.g. Fj, see §5), and 
the behavior of rjnfl is set by the disk temperature in the 
outermost, directly irradiated and often optically thin 
regions. This motivates more careful treatments of disk 
thermodynamics in the future. However, we believe that 
this simplification may affect our results at most at the 
level of tens of per cent. 

We expect our results on the global disk evolution to be 
rather insensitive to the exact form of the binary torque 
density A(r). However, equation (16) shows that heating 
of the innermost disk regions is directly determined by 
A(r), for which we have adopted a simple form (3). Fu¬ 
ture treatments of the circumbinary disk physics would 
benefit from the more accurate torque density descrip¬ 
tion (Goodman & Rafikov 2001; Rafikov 2002; Rafikov 
& Petrovich 2012; Petrovich & Rafikov 2012). 

Our models have typically been evolved for 5 Myr with¬ 
out any external mass loss from the system. In practice, 
disk is likely to be losing mass due to photoevaporation 
(Alexander 2012), which will modify its structure and 
evolution at late times. 

8.7. Comparison with other work. 

Several authors have previously addressed different as¬ 
pects of the circumbinary disks and we discuss how our 
results fit in the context of these studies. 

Alexander (2012) and Martin et al. (2013) presented 
evolutionary circumbinary disk models based on equation 
(1), focusing on the role of photoevaporation and dead 
zone, respectively. Both studies have chosen the torque 
density prescription in the form (3) with / = 0.5. This 
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value is much larger than what we use in our work, result¬ 
ing in wider inner cavities in these studies — (4 — 5)ab, 
which is « 2 times larger than suggested by the results 
of hydrodynamical simulations (MacFadyen & Milosavl- 
jevic 2008; Pelupessy & Portegies Zwart 2013). Because 
our disks extend closer to the central binary (down to 
~ 2 ab), tidal and viscous heating drive inner disk tem¬ 
perature to significantly higher values in our models, up 
to 2000 K as opposed to 800 K in non-accreting models 
models of Martin et al. (2013) with the same disk mass. 

We also provide a more accurate treatment of the 
disk thermodynamics. For comparison, Alexander (2012) 
have simply prescribed a fixed temperature profile in the 
disk. Martin et al. (2013) provide a treatment of disk 
thermodynamics similar to ours and account for viscous 
and tidal heating, albeit with a simplified opacity behav¬ 
ior K oc ® adopted for the whole disk. However, their 
prescription for irradiation heating assumes Tirr cx; 
which is not appropriate for the externally irradiated, 
flared disk. As a result, the outer regions of their disk 
models receive much less heating than the flared irradi¬ 
ated disks, which should have implications for the viscous 
spreading and global evolution of their disks. 

Clanton (2013) explored the issue of iceline location in 
a circumbinary disk and found rice to lie rather close to 
the binary. However, this calculation did not account for 
the tidal dissipation, which as we show in this work, dom¬ 
inates heating in the inner disk. Shadmehri & Khajenabi 
(2015) have improved the rice calculation by accounting 
for tidal heating, and found significantly larger values 
of Tice. However, both of these studies did not prop¬ 
erly account for the changes of the disk structure caused 
by the injection of the angular momentum at its center. 
Instead, they simply adopted E(r) structure that would 
correspond to a standard constant-M disk with no torque 
at the center. As we show in our work, this assumption 
is not justified in general and S(r) and Fj{r) distribu¬ 
tions should be quite different in realistic circumbinary 
disks from their analogues in conventional constant-M 
circumstellar counterparts. 

9. IMPLICATIONS FOR PLANET FORMATION 

We now discuss the ramifications of our results for cir¬ 
cumbinary planet formation. 

Suppression of accretion (even partial) by the binary 
torque implies that the circumbinary disk retains more 
mass that its circumstellar counterpart at the same mo¬ 
ment of time. In fact, there is some observational evi¬ 
dence in support of this conjecture. Harris et al. (2012) 
hnd that circumbinary disks tend to be significantly 
brighter at sub-mm wavelength than their counterparts 
around the single stars, although one has to keep in mind 
that the size of their sample of circumbinary disks with 
detected emission is small (only four systems). Higher 
disk mass persisting for longer has a number of impor¬ 
tant consequences for planet formation in circumbinary 
disks. 

First, higher surface density of gas and solids should 
speed up the growth of planetesimals due to larger mass 
supply. For illustration, comparison of Figures 6 and 9 
demonstrates that e.g. at 3 AU E in a circumbinary disk 
is higher by ~ 3 compared to the circumstellar disk of 
the same initial mass M 4 = 0.05Mc at t = 0.5 Myr. This 


implies that in the former case planetesimal growth by 
coagulation would be faster by roughly the same factor. 

Second, gravitational effect of a more massive disk sig¬ 
nificantly reduces eccentricity excitation by the binary 
gravity (Rafikov 2013a; Silsbee & Rafikov 2015). As 
shown by (Silsbee & Rafikov 2015), this lowers rela¬ 
tive speeds with which planetesimals collide, alleviating 
the so-called fragmentation problem for planet formation 
in binaries (Meschiari 2012; Paardekooper et al. 2012; 
Marzari et al. 2013). 

Third, the isolation mass of the cores forming in mas¬ 
sive circumbinary disks should be larger than in their 
circumstellar counterparts, since it scales as This 

may allow direct formation of the planetary cores mas¬ 
sive enough to trigger runaway gas accretion and turn 
into giant planets. 

Our results have non-trivial implications for another 
issue relevant for reducing collisional velocities of plan¬ 
etesimals in binaries, namely the precession of the cen¬ 
tral binary. It was shown by Rafikov (2013a) and Sils¬ 
bee & Rafikov (2015) that disk gravity often dominates 
precession of the eccentric central binary, which in turn 
reduces planetesimal eccentricity excitation by the non- 
axisymmetric component of the binary gravity. The rate 
of binary precession is set mainly by the disk surface den¬ 
sity at its inner edge, which was previously estimated 
using simple models of passively heated (externally irra¬ 
diated) disks (Rafikov 2013a; Silsbee & Rafikov 2015). 
However, our results show (see Figures 6-8) that because 
of the tidal and viscous heating T is considerably higher 
in the inner disk regions than equation (22) would pre¬ 
dict, implying higher viscosity and lower E at the cavity 
edge compared to the prediction (25). This may lower 
the precession rate of the central binary by a factor of 
several compared to the existing models. 

On the other hand, outward displacement of the iceline 
by the vigorous tidal and viscous heating may negatively 
affect core assembly in the inner (within few AU) regions 
of circumbinary disks by reducing the surface density of 
solids (per E of gas) there. This effect may favor forma¬ 
tion of the observed Kepler circumbinary planets in the 
outer regions of their parent disks (e.g. outside ~ 5 AU) 
with subsequent inward migration into their present day 
location, a scenario that is preferred on other grounds as 
well (Paardekooper et al. 2012; Rafikov 2013a; Silsbee & 
Rafikov 2015). 

10. SUMMARY. 

We explored properties and viscous evolution of pro¬ 
toplanetary disks around stellar binaries, which are the 
birthplaces of circumbinary planets such as the ones de¬ 
tected by the Kepler mission. Our main conclusions are 
briefly summarized as follows. 

• Circumstellar disks are in many ways different from 
their circumstellar analogues (§6.3). Suppression 
of accretion onto the central object by the binary 
torques not only allows the disk to retain mass in 
the disk for longer. Tidal barrier imposed by the 
central binary modifies the basic character of the 
radial distribution of the viscous angular momen¬ 
tum flux Fj{r), leading to a flat Fj profile in the 
case of no accretion (§3.2). This has important 
implications for the disk structure — circumbinary 


disks contain more mass and release more energy by 
viscous heating in their central regions than their 
circumstellar counterparts of the same mass. 
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• Very importantly, we derive a set analytical rela¬ 
tionships for the viscous evolution of disk proper¬ 
ties (§5) that are verified and calibrated by the de¬ 
tailed simulations with realistic inputs (§6). These 
relations form a basis for quantitative understand¬ 
ing of the role of different parameters of the system 
(disk mass, viscosity, etc.) on its evolution. They 
are flexible and can be generalized to account for 
additional effects such as the non-negligible accre¬ 
tion onto the central binary (§8.2). 

• Dissipation of the binary-driven density waves 
dominates heating of the inner disk, within 1-2 
AU (§8.1). This energy source (absent in disks 
around single stars) raises inner disk temperature, 
and pushes the iceline further out (to ^ (5 — 10) 
AU) compared to circumstellar disks (§8.5). Ir¬ 
radiation by central binary starts to control disk 
temperature only outside ^ 8 AU. 

• SED of a circumbinary disk is different from SED 
of its circumstellar counterpart of the same mass. 
Dissipation of the density waves gives rise to a dis¬ 
tinctive bump in the SED around lO^m that may 
facilitate identification of circumbinary disks when 
the binarity of the central source is not obvious 


(§7). 

• Viscous and tidal energy release in the central re¬ 
gion give rise to self-shadowing of the disk by its 
inner parts out to ~ 20 AU (§7.1). 

• Tidal coupling to the disk continuously removes an¬ 
gular momentum from the central binary, shrinking 
its orbit and potentially resulting in its merger into 
a single star (§8.3). 

• Circumbinary disks are in many ways more favor¬ 
able sites of planet formation than their analogues 
around single stars (§9). This is in agreement the 
occurrence rates of circumbinary planets inferred 
from statistics of Kepler systems (Armstrong et al. 
2014). 

Our study thus provides a basis for understanding the 
systems in which Kepler circumbinary planets were born. 
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APPENDIX 


NUMERICAL SETUP FOR CIRCUMBINARY DISK EVOLUTION. 

To explore viscous evolution of a protoplanetary disk we numerically solve equation (1). We first recast it as a mass 
conservation law with an advection term, and then use the finite-volume method to discretize it with an explicit-upwind 
scheme. The time step size is determined adaptively by computing the minimum time for any cell in the grid with 
nonzero mass to be emptied if the fluxes are held constant; then a small fraction (0.01%) of that time is used as the 
time step size. We use the software package FiPy (Guyer et al. 2009) for the discretization and the integration of 
equation (1). Our source code is freely available®. 

We use a logarithmically spaced radial grid that spans from =0.2 AU to Vgut = 2 x 10^ AU with 150 grid points 
(we verified that our results are converged at this resolution). We consider a central binary with g = 1, Me = Mq 
(i.e. Mp = Ms = O. 5 M 0 ), flb = 0.2 AU, Le = 2Lq, and three disk masses Md = 0.01 Me, 0.05 Me, and 0.1 Me- The 
viscous a-parameter is set to a = 0.01 in all our runs. 

Disk thermodynamics enters equation (1) only via the a-parametrization of viscosity (2), and is self-consistently 
treated by solving equation (19) at each time step. Our temperature prescription thus allows for viscous, tidal, and 
irradiation heating and interpolates over optically-thick and -thin limits. 

Disk irradiation by the binary is treated in a simplified way. Namely, we assume that the incidence angle ^ is always 
that of a purely irradiation-dominated disk. For that reason is not computed using equation (15) with the actual 
T{r) setting h{r). Instead, C is simply given by equation (24). As the disk may become non-flared at certain radii, 
shadowing exterior portions of itself (which we indeed find to be the case, see §7.1), the irradiation may be absent 
in some parts of the disk. For simplicity, we neglect this self-shadowing effect and discuss the ramifications of this 
approximation in §7.1. 

To estimate optical depth, we adopt the opacity fits from Zhu et al. (2009), which assume k behavior in the form 
K = KiT°‘'p^' in a number of regions in T-p space. Parameters Ki, a^, bi depend on the physical condition in the i-th 
opacity regime, and change as dust particles of different composition sublimate, molecules are dissociated, and different 
species are ionized (see Zhu et al. (2009) for details). This k prescription spans from electron scattering opacity at 
high T to opacity dominated by water ice grains at low T. 

As our initial condition we start the disk as a narrow ring of mass at some initial radius tq and let it evolve viscously. 
The initial radial distribution of E is given by a Gaussian ring of width ar rp 


E(r-,t = 0) 


Md 

(27r)®/2ro(Tr- 


{r - rpf 
2 a e 


(Al) 


where we took ro = 20 AU and ar = 2 AU. After several viscous timescales our results should not depend on this 
initial condition. 

We use the following boundary conditions (BCs). At the outer boundary of our computation domain we set mass 
accretion rate to zero: 


dFj 

dl 



M{rout,t) = 0 . 


(A2) 


In practice, our spreading disks never reach the outer boundary so that this BC is well observed. At the inner boundary, 
in a circumstellar case, we impose a uniform (in radius) accretion rate BC via 


dFj 

dl 



Fj{rin) 

7. 

('in 


(A3) 


where Im = l{rin) is the specific angular momentum at the inner boundary. This form of BC is motivated by equation 
(11). Note that we do not specify the actual value of M, instead, we allow M{rin) to be self-consistently determined 
by the viscous stresses. In a circumbinary case, because of high masses of the stellar components, A in the form (3) 
results in vanishingly small E(ri„) and M(rm), resulting in Mb = 0 (see §3.1), so that the inner BC is irrelevant. In 
real disk, M{rin) may be different from zero because of non-axisymmetric effects not captured by our ID approach 
(D’Orazio et al. 2013). For simplicity we do not consider this possibility here. 

The injection of angular momentum in the disk center by the binary is modeled as follows. For |r — ab\ > h{r), we 
follow the prescription (3), while for \r — ab\ < h{r) we use A(r) oc r, continuously matching at r = ± h(r). Note 
that this introduces a discontinuity in the derivative of A(r) at this point. To avoid numerical problems we remove 
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the discontinuity by smoothing A. Since we are using a high mass ratio {q = 1) and matter never gets close to the 
binary, the details of the smoothing are unimportant. 

Moreover, as previously stated in §3, the details of the tidal disk-binary coupling occurring at r < ta are irrelevant 
for the global disk evolution at r > ta, as long as the central angular momentum source gives rise to a proper inner 
BC. This is how we exploit the prescription (3), not worrying about the details of A(r) but simply making sure that 
our ID disk profile reproduces gross features of the more detailed 2D numerical calculations of circumbinary disks. As 
a particular metric for comparison we use the size of the inner cavity, which MacFadyen & Milosavljevic (2008) found 
to be about 2ab. In all our calculations we tune the value of the parameter / in the expression (3) so that in each 
of our runs cavity edge is at rj, ss 2af,. This typically results in / « 10“^-2 x 10“^, with specific values for each disk 
model indicated in Figures 6-8. 


SPECTRAL ENERGY CALCULATION 


The energy flux J-, emitted by a unit surface area element of the disk at each radius r, consists of several contributions. 
First, there is a (one-sided) flux due to viscous and tidal dissipation Ty and J);id given by equations (13) and (16). 
Second, part of the incoming stellar radiation is intercepted by the superheated dust layer, re-emitted towards the disk, 
absorbed and then lost to space. Finally, the remaining part of stellar radiation intercepted by superheated grains is 
re-emitted directly to space. The relative partition between the last two contribution depends on the optical depth of 
the disk. 

The first two flux contributions are re-emitted by the disk at the characteristic temperature Te, which we compute 
according to the following formula: 

aT^ = (1 + T-l) (Ty + J-tid) + Tirr- (Bl) 

In the optically thin limit this expression reduces to the midplane disk temperature given by equation (17). 
Superheated dust layer emits at the temperature Tgh given by (Chiang & Goldreich 1997) 


Tsb{r) 


Lc 


1/4 


167ro-e(Tsh)?’2 


(B2) 


Here e(Tsh) = Qabs(7sh)/Qabs(7*) is the ratio of absorption efficiencies of small grains at Tsh and stellar temperature 
r*, for which we adopt the approximation of Chiang & Goldreich (1997) e{T) T/Ti,. 

The one-sided disk SED is then computed using the following interpolating formula, which is designed to correctly 
reproduce the limiting cases of optically thick and thin disk: 

F, = + (f^) (B3) 

For simplicity we do not consider t to be a function of frequency in this expression. When integrated over if this 
expression results in the total flux (J-y + J^tid) + for any t, in agreement with energy conservation (recall that 
Jdrr is a half of the stellar flux incident on a disk surface). Note that J^rr/ = 2e(Tsh)C) see equations (14) and 

As stated in §4.2, we are ignoring the difference in values of opacity for the radiation at temperatures Tg and Tsh. 
In the optically thick limit (r ^ 1) the amount of flux J^y + T’tid + -Tirr is emitted by the disk at effective temperature 
Te, while another Tirr is radiated by the superheated layer at the temperature Tsh- In the limit of small optical depth 
(r <C 1) disk emits only the sum of viscous and tidal energy fluxes Ty -I- Tfid, while the two superheated layers (both 
of which are visible because the disk is transparent as r —> 0) radiate 2 J-iyy. 

Note that, as described in Appendix A, in computing the SED we assume the grazing incidence angle of the stellar 
radiation to be always given by equation (24). This assumption may not be justified in intermediate radii, which can 
be shadowed, as discussed in §7.1. Also, use of Tlrr in the form (14) at all distances neglects the fact that superheated 
dust layer becomes inefficient at absorbing stellar radiation at large radii. Neglect of this issue artificially boosts the 
emission produced by the superheated layer compared to the disk emission, see §7. 








